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FOREWORD 


A  multiple  zone  computational  method  applicable  to  finned  bodies  in 
supersonic  flight  is  described.  At  each  cross-section,  the  computational 
domain  is  broken  into  a  number  of  zones,  each  of  which  can  be  described  using 
a  simple  mapping.  The  zones  are  chosen  to  insure  that  the  bow  shock,  fin  and 
body  surfaces  coincide  with  zone  edges.  At  zone  interfaces,  the  mesh  may  be 
discontinuous  in  a  direction  normal  to  the  edge.  Along  zone  edges surfaces 
are  allowed  to  form  or  disappear  during  the  computation.  The  surface  edges 
which  form  at  points  where  surfaces  appear  or  disappear  are  treated  using  a 
local  analysis.  The  resulting  algorithm  is  applicable  to  a  wide  range  of 
configurations  including  interior  as  well  as  exterior  supersonic  flow  fields. 
Calculations  are  compared  to  experiment  for  several  different  missile  shapes. 
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1.  INTRODUCTION 

A  practical  means  of  predicting  the  inviscid  flow  field  about  a 
supersonic  missile  is  to  numerically  solve  the  steady,  three  dimensional  Euler 
equations.  As  is  shown  in  Figure  1,  such  an  approach  is  started  at  a  cross- 
sectional  plane  near  the  missile  nose  tip,  using  a  known  flow  field  which  can 
be  obtained  either  from  a  conical  or  blunt  body  solution.  The  hyperbolic 
nature  of  Euler  equations  in  supersonic  flow  allows  this  flow  field  to  be 
marched  in  the  axial  direction  as  long  as  the  flow  remains  supersonic 
everywhere. 

A  major  problem  encountered  in  steady,  supersonic  calculations  is  the 
treatment  of  complex  cross-sectional  geometries  which  arise  in  conjunction 
with  configurations  featuring  fins,  inlets  and  tails.  Most  inviscid 
treatments*"'  for  arbitrary  bodies  employ  a  wrap-around  transformation  which 
maps  both  the  body  and  fin  surfaces  onto  the  edge  of  the  computational  regions 
as  is  shown  in  Figure  2.  Mapping  methods  have  been  developed  for  relatively 
arbitrary  configurations,-  including  bodies  with  wings.  However,  for  missiles 
with  multiple  wings  and  fins  which  are  thin  and  have  sharp  edges,  this  type  of 
approach  has  several  drawbacks.  Solutions  may  be  sensitive  to  small  changes 
in  the  transformation  and  it  is  often  difficult  to  control  mesh  spacing 
throughout  the  flow  field.  Also,  a  significant  portion  of  the  computation  may 
be  associated  with  constructing  the  computational  mesh. 

An  alternative  approach  for  treating  complex  missile  and  aircraft 
configurations  is  a  multiple  zone  strategy.  Here,  the  computational  domain  in 
each  cross-section  is  divided  into  several  zones,  each  of  which  can  be  treated 
using  a  simple  mapping.  A  sample  zone  structure  for  a  body-wing-tail 
configuration  is  shown  in  Figure  3.  This  avoids  many  of  the  mapping 
complexities  associated  with  wrap-around  meshes;  however,  it  introduces 
additional  program  complexity  associated  with  interconnecting  adjacent 
zones.  Despite  these  difficulties,  the  described  multiple  zone  method  is 
capable  of  treating  a  wider  variety  of  configurations  than  a  wrap-around  mesh 
method.  Furthermore,  use  of  a  separate  mesh  in  each  zone  allows  selected 
portions  of  the  flow  field  to  be  resolved  in  much  greater  detail. 

Multiple  zone  methods  have  been  used  in  applications  of  the  unsteady 
Euler  equations  to  a  number  of  complex  configurations  including  inlets  with 
blunt  lips,  airfoils  with  flaps,  and  body-wing  combinations.  1 0  In  these 
studies,  two  basic  multiple  zone  structures  have  been  used;  overlapping  zones, 
and  abutting  zone,  which  only  interface  along  a  common  boundary.  Regardless 
of  the  type  of  zone  interface  employed,  information  is  transmitted  between 
zones  by  interpolation.  If  shocks  are  to  be  captured,  this  interpolation 
procedure  should  preserve  the  global  conservation  property  of  the  scheme. 


1 


NSWC  TR  85-484 


This  paper  focuses  on  the  application  of  a  multiple  zone  approach  to 
missiles  and  aircraft  with  thin,  sharp  fins  or  wings.  Such  shapes  are 
conveniently  treated  using  abutting  zones  with  portions  of  the  boundaries 
taken  to  coincide  with  segments  of  the  body,  fin  or  tail  surface,  and  the  bow 
shock.  Since  the  cross-sectional  geometry  and  bow  shock  location  vary  along 
the  missile  axis,  the  transformation  in  each  zone  must  be  recomputed  at  every 
step  in  the  calculation.  To  promote  efficiency,  simple  analytic  mappings  are 
used.  The  complexity  of  the  interfacing  algorithm  is  reduced  by  requiring 
that  points  along  adjacent  zone  boundaries  be  coincident.  This  allows  a 
conservative  edge  point  treatment  to  be  applied  which  does  not  require 
interpolation.  The  mesh  may  be  discontinuous  normal  to  the  zone  edge, 
however,  this  may  reduce  the  local  solution  accuracy.  Each  zone  is  integrated 
using  the  explicit  MacCormack  scheme  with  wall  and  shock  points  being  treated 
by  characteri sti c  analysis.  Special  computational  methods  are  used  at 
corners,  body-fin  junctures,  surface  edges  and  surface  slope  discontinuities. 

A  Multi-zone,  Supersonic  Euler  code  (MUSE)  has  been  developed  which 
implements  the  multiple  zone  procedure.  The  computational  approach  is  similar 
to  that  of  the  SWINT  code  reported  on  in  References  12-14.  However,  Euler 
equations  are  solved  in  cartesian  rather  than  cylindrical  coordinates.  This 
simplifies  the  analysis  of  points  lying  on  the  zone  edges  and  results  in  a 
system  of  equations  which  is  in  strong  conservation  form.  In  principal,  this 
method  can  be  applied  to  any  internal  or  external  supersonic  flow  whose  cross- 
section  can  be  described  by  a  number  of  quadrilateral  zones.  The  MUSE  code 
can  be  applied  to  more  general  configurations  than  the  SWINT  code  and  can 
account  for  fin  thickness.  However,  the  MUSE  code  also  requires  30%  to  60% 
longer  to  run  per  mesh  point,  per  time  step  and  uses  a  similar  increase  in 
storage  for  fine  mesh  runs.  On  configurations  which  can  be  handled  by  both 
SWINT  and  MUSE,  similar  results  are  obtained. 

2.  MULTIPLE  ZONE  METHODOLOGY 

The  multiple  zone  procedure  is  implemented  by  dividing  the  shock  layer 
into  non-overlapping,  quadrilateral  zones.  Zone  edges  are  taken  to  coincide 
with  walls  and  the  bow  shock,  if  one  exists.  Here,  a  wall  refers  to  any 
surface  at  which  a  tangent  flow  boundary  condition  is  applied,  such  as  body, 
fin  or  tail  surfaces.  Points  interior  to  the  shock  layer  may  also  lie  on  zone 
edges.  Juncture  points,  which  are  the  intersection  of  two  wall  surfaces,  must 
occur  at  zone  corners.  Points  lying  on  zone  edges  may  change  their  type  at 
different  axial  stations  from  interior,  to  wall,  and  vice  versa.  A  similar 
transition  between  interior,  wall  and  juncture  points  can  also  occur  at  a  zone 
corner.  This  allows  the  simulation  of  inlet  lips  and  wing  leading  or  trailing 
edges.  In  the  present  procedure,  the  intent  is  to  fit  the  bow  shock  to 
provide  a  boundary  for  the  computational  domain  and  to  let  the  conservative 
properties  of  the  different  scheme  capture  imbedded  shocks.  Figure  3 
illustrates  an  example  of  a  four  zone  model  at  a  cross-section  featuring  a 
body,  wing  and  tail. 

All  interior  points  are  governed  by  the  steady  Euler  equations,  which  are 
cast  in  the  cartesian  coordinates  (x,y,z)  defined  in  Figure  1.  In  each  zone, 
these  equations  are  transformed  to  computational  space  (c,n,t)  using  the 
procedure  outlined  below.  The  transformed  equations  are  in  conservation  form 
and  are  written  as: 
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3U  +  9_F  +  3G 

3?  35  3n 

where 


0 


U  =  (U/J);  F  =  (5zU  +  5XF  +  CyG)/J ; 

G  =  (nzU  +  nxF  +  nyG)/J;  J  =  5xny  -  5ynx  (1) 


“pw  ' 

”pu 

■pv 

p+pw 

;  F  = 

puw  2 

;  G  = 

pVW 

puw 

p+pu 

pvu  2 

_pVW 

pvu 

_p  +  pv 

and  p,  p,  u,  v,  w  are  the  pressure,  density  and  cartesian  velocity  components, 
respectively.  These  equations  are  closed  using  the  constant  total  entropy 
constraint  and  the  perfect  gas  equation  of  state: 


P  ”Y 

/ --T +  4  (u2+v2+w2)  =  H  (constant)  =  7 — —  +  4  (u^  +  v2  +  w2) 
(y-l)p  2  o  (y-l)pOT  2  °° 

A  different  set  of  computational  coordinates  is  used  in  each  zone. 


(2) 


a .  Zone  Description 

In  the  crossflow  planes  z  =  constant,  zones  are  generalized 
quadrilaterals  as  shown  in  Figure  4.  Each  zone  is  described  with  respect  to 
( s ,  x ,  v )  coordinates,  which  either  represent  cylindrical  (r,c}>,z)  or  cartesian 
(x,y,z)  coordinates.  Cylindrical  coordinates  facilitate  treatment  of  wing- 
body  configurations,  while  cartesian  coordinates  can  be  applied  to  wing  alone 
cases.  The  numbering  system  used  to  designate  the  edges  and  corners  of  a  zone 
are  indicated  in  Figure  4.  The  locations  of  edges  1  and  3  are  defined  by  the 
functions  b(x,v)  and  c(x,v),  while  the  coordinates  of  edges  2  and  4  are 
defined  by  ^(s.v)  and  o{ s.v),  respectively.  When  (s,x,v)  represent 
cylindrical  coordinates,  the  bow  shock  may  be  fitted,  but  it  must  be  located 
on  edge  3. 


The  mesh  in  each  zone  is  defined  by  the  transformation 
T-l^,  where  T^ :  (5,n,?)  (s,x,v)  and  T 2  (s,x,v)  (x,y,z).  T2  is  given  by: 


x  =  s  cos  (x) 
y  =  s  sin  (x) 
z  =  v 

x  =  s 
y  =  x 
z  =  v 


cylindrical  coordinates 

(3) 

cartesian  coordinates 


while  T^  is  defined  as. 


5  =  b(x' ,?)  +  [c(x",c) 
t  =  ct(s‘  ,?)  +  0(s" ,?) 
v  =  ? 


b(r*,?)]fU) 
o-(s'  5)]g(n) 


(4) 
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where: 

T  '  =  T 4(e  )  +  (t  ]  (c  )  -  t4(c  )  )g(n  ) 

T"  =  T^c)  +  (xh?)  -  T^(c))g(n) 

S'  =  sffo  )  +  (s~(c  )  -  s^(c))f(c) 

S"  =  sj(c)  +  (s|(c)  -  sj(c))f(5) 

and  (s-,x-),  i  =  1,2, 3, 4  are  the  coordinates  of  the  corners.  Here  f  and  g  are 
mesh  clustering  functions  with  f(0)  =  g(0)  =  0  and  f(l)  =  g(l)  =  1.  The 
computational  domain  for  each  zone  is  bounded  by  0  <  g  <  1,  0  <  n  <  1,  and 

g  >  0 . 


The  metric  quantities  E,5,£,n»n»n»  appearing  in  Equation  (1) 
must  be  evaluated  at  every  point^in  the  calculation.  In  addition,  at  wall 
points  on  edges  1  or  3,  the  quantities  g  x  ,  g  ,  g  ,  are  required. 

Similarly,  n  ,  n  ,  n  ,  must  be  calculated  at  wan  points  on  edges  2  and 
4.  Analytic^xpretsio^s  for  these  tenns  are  provided  in  Appendix  A. 


Alternatively,  g  ,  g  ,  g  ,  n  ,  n  ,  n  can  be  evaluated  numerically.  This 
is  accomplished  by  determining  r  *  g  Z,  g  using  central  differences  and 
then  computing  the  metric  quantiti es^ f roffi :  ^ 


5x  =  <\V  v«,/j;  nx  =  (-Ts5x  +  Vx>/J 
{y  *  (Vy  '  Vy)/j;  ny  =  ('Vy  +  Vy)/J 

?z  =  (snT?  ”  Tr)Sc^’  nz  =  “  S£x,r)/J 


where : 


J  =  Vn  -  Vc 


Here  and  s^  are  evaluated  analytically, 
b.  Zone  Linking 


(5) 


Information  must  be  passed  between  zones  connected  by  abutting  edges. 

Two  edges  are  abutting  if  they  possess  segments  which  are  coincident  in 
cartesian  space  and  not  separated  by  a  zero  thickness  surface.  To  simplify 
the  interfacing  scheme,  the  zones  are  defined  to  ensure  that  abutting  edges 
satisfy  the  following  rules: 

1.  The  two  sets  of  points  on  each  of  the  abutting  edge  segments  must  be 
coincident  in  cartesian  space. 

2.  Odd  numbered  edges  must  abut  to  odd  numbered  edges  and  even  to  even. 
Also,  edges  with  the  same  number  cannot  abut. 

3.  A  zone  corner  which  is  on  an  abutting  edge  segment  can  only  be 
coincident  in  cartesian  space  with  other  zone  corners. 

Rule  one  requires  that  meshes  in  zones  which  abut,  produce  coincident  grid 
points  in  cartesian  space  along  abutting  edge  segments.  Restrictions  are  not 
imposed  on  the  mesh  spacing  normal  to  abutting  edges  and  abrupt  changes  in 
mesh  size  can  occur. 
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The  computational  algorithm  treats  each  zone  separately,  except  along 
abutting  edges.  Here, 'property  and  metric  information  from  the  abutting  zones 
must  be  passed.  Generation  of  the  appropriate  metrics  is  accomplished  by 
viewing  the  mesh  in  the  abutting  zones  as  a  single,  extended  computational 
space.  For  example,  consider  the  case  of  two  abutting  zones,  I  and  II,  with 
computational  spaces  of  (5  1  ,n  '  ,4  1 )  and  U",n"»s")  respectively.  Apply  the 
discontinuous  transformation: 


s  =  U 1 

in 

zone 

i;  k 

"  + 

Cq  in  zone 

n} 

n  =  {n  1 

i  n 

zone 

i;  Bn 

"  + 

in  zone 

n} 

4  =  {?' 

in 

zone 

i;  4" 

in 

zone  11} 

where  a  =  A£;,/A£",  B  =An'/An".  The  constants  Cq  and  c^  can  be  chosen  so  that 
the  two  zones  correspond  in  (g,n,c)  to  a  square  and  a  rectangle  joined  along 
the  abutting  edge.  Rule  1  implies  that  the  combined  mesh  is  continuous  in 
cartesian  space  along  the  abutting  edge  and  has  a  uniform  spacing  A£  ,An  in 
computational  spaces  (£  ,n  ,4  ) .  Although  analytic  values  of  x  and  v  are 
continuous  along  the  abutting  edge,  those  of  x  ,x  ,4  ,4  may^be  discontinuous. 
These  discontinuities  are  removed  in  the  present  $aper  By  defining  single 
values  for  these  derivatives  using  standard  centered  differences  in  5  ,n 
variables.  Then  along  the  abutting  edges, 


4 


5 


x 


4  '  =  4  "/B  =  4 
n  n  n 


x  1  =  x  "/b  =  X  , 

n  n  n 


(6) 


where  the  overbar  denotes  values  from  centered  differences.  The  required  metric 
quantities,  %  ,  5  ,  £  ,  n  ,  n  ,  nz,  are  then  calculated  using  Equation  (3). 

As  will  be  shown  in  tne  next  -section,  the  above  modifications  to  the  metrics 
at  abutting  edges  permits  the  application  of  the  MacCormack  scheme  with  only 
slight  modification.'  Interpolation  is  not  required  and  the  scheme  satisfies  a 
global  conservation  property. 


The  use  of  central  differences  to  calculate  metrics  along  abutting  edges 
is  consistent  with  a  second  order  accurate  scheme  when  the  mesh  variation 
across  the  abutting  edge  is  smooth  in  cartesian  space.  In  other  cases,  a 
degeneration  in  accuracy  can  be  expected.  In  Figure  5,  results  obtained  for 
flow  over  a  cone  is  shown  using  a  two  zone  model  which  features  a  5:1  jump  in 
mesh  spacing  across  the  abutting  edge.  The  abutting  zone  edges  are  located  in 
the  region  of  smooth  flow  and  good  agreement  is  obtained  with  a  uniform  mesh 
calculation.  When  the  abutting  edge  is  in  a  region  of  flow  featuring  a  shock, 
the  shock  trajectory  is  accurately  captured  in  both  zones,  but  additional 
noise  is  introduced  on  the  downstream  side  of  the  shock.  Figure  6  illustrates 
such  a  case  using  a  biconic  configuration. 


3.  NUMERICAL  PROCEDURE 


The  solution  in  each  zone  is  advanced  in  its  own  computational  space.  A 
common  marching  step  size  is  employed  in  all  zones,  which  is  determined  by 
applying  the  CFL  condition  to  all  points  in  the  calculation. 
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a.  Interior  Points 

Interior  points  are  advanced  using  the  explicit  MacConnack  scheme  applied 
to  Equation  (1).  This  is  a  predictor-corrector  scheme  which  can  be  written  as 
fol lows: 


V 

j 

!  ui.jn  -  (Fi+k.j  - F 

)  ££  - 
i+k-1 ,j ;  A5 

<Gi 

j+* 

Gi ,  j+a - 1 ' 

A_£ 

An 

ui,j - 

1  Cui 

,j 

+  Ui  J  ”  ^i-k+1  ,j 

,  j -£  +  l 

A? 

An 

where: 

k  = 

1 

forwards-backwards 

di fferenci ng 

i  n 

the  5 

di  recti  on 

0 

backwards- forwards 

di fferenci ng 

in 

the  5 

di  recti  on 

i  = 

1 

forwards-backwards 

di  fferenci  ng 

in 

the  n 

di  recti  on 

0 

backwards-forwards 

di fferencing 

in 

the  n 

di recti  on 

(7) 


Here  the  bar  superscript  indicates  fluxes  evaluated  using  predictor  values. 
The  same  varient  of  this  scheme  (i.e.,  forwards-backwards  or  backwards- 
forwards)  is  applied  to  each  zone. 

The  step  size  A E.  is  chosen  using  the  CFL  condition,  which  requires  that 
the  following  condition  be  satisfied  at  every  point: 


A? 


2  2 
(w  -a  )An 


maxJy^  ,y  2>u 


(8) 


where: 

vl  =  M  |wA-a2?zj  +  a  V(A  -  w?z)2  +  (w2  -  a2)[sx2  +  ?y2] 
u2  =  | wB  -  a2nz|  +  a -f/(B  -  wnz)2  +  (w2  -  a2)  [nx2  +  ny2] 
u 3  =  |wB  -  a2nz  +  (-l)j6(WA  -  a25z)  | 

+  a  ^/[B-wn  z+(-l)^6 (A-w?z)j2+(w2-a2)[(nx+(-l)J65x)2+(ny+(-l)J6?y)2] 


6  =  An /A? 

j  =  M0Dj_k+£  ,2) 

A  =  V?  q 
B  =  Vn  q 
q  =  (u,v,w) 

At  the  completion  of  every  step,  it  is  necessary  to  decode  the 
conservation  quantities  to  determine  p,p,u,v,w.  This  is  accomplished  as 
fol lows : 
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u  =  U3/U1 

v  =  U4/U1  (9) 


./II  _  _ 

W  =  TT+lT  1  +  y\  -  [2Hq 

L  Y 

ui 

U2 

2 

U3 

U2 

2 

U4 

U2 

- 1 

f - 1 

CM 

p  =  (U2  -  wllj  )J 
p  =  U^J/w 


At  interior  points  lying  on  an  abutting  edge,  the  numerical  scheme  must 
be  modified  to  allow  access  to  information  in  the  abutting  zone.  This  is 
accomplished  by  viewing  the  two  abutting  zones  as  part  of  an  extended  space, 
as  was  done  in  Section  2.  In  this  space,  Equation  (1)  holds  both  in  and  on 
the  edges  of  the  abutting  zones.  The  metric  treatment  of  Section  2  provides  a 
continuous  definition  of  metric  quantities  and  the  fluxes  U,F,G  in  the 
neighborhood  of  the  abutting  edges.  This  allows  the  MacCormack  scheme  to  be 
applied  to  Equation  (2)  at  points  located  along  the  abutting  edge.  To 
accomplish  this,  it  is  necessary  to  convert  the  fluxes  to  the  extended  space 
coordinates.  For  example,  to  advance  a  point  in  zone  I  which  abuts  zone  II, 
the  quantities  appearing  in  the  MacCormack  scheme  are  defined  as  follows: 


A5  =  J 

A?  ' 

i  n 

zone 

1, 

^  11 

dAE, 

i  n 

zone 

11 

An  =  • 

An  ' 

i  n 

zone 

1, 

aAn" 

in 

zone 

11 

U  =  ■ 

U' 

i  n 

zone 

1, 

U"/ab 

in 

zone 

ii 

G  = 

G' 

i  n 

zone 

1, 

G7a 

in 

zone 

11' 

F  = 

'  F' 

i  n 

zone 

1, 

F"/b 

in 

zone 

ir 

Here  '  and  "  denote  quantities  evaluated  in  zone  I  and  zone  II  computational 
spaces  respectively.  An  analogous  procedure  is  used  to  advance  an  edge  point 
in  zone  II.  The  above  methodology  for  edge  points  ensures  that  the  decoded 
flow  variables  at  abutting  points  in  cartesian  space  are  identical  at  the  end 
of  both  the  predictor  and  corrector  steps.  Interior  points  lying  at  zone 
corners  are  treated  in  an  analagous  fashion.  Flere  the  extended  space  spans 
the  four  zones  adjacent  to  the  corner. 

b.  Mall  Points 

The  governing  equations  are  cast  in  characteristic  form  and  differenced 
in  the  manner  proposed  by  Kentzer.^  One  of  these  characteristic  equations  is 
inadmissable  since  it  carries  information  from  outside  the  computational 
domain,  and  it  is  replaced  by  the  tangent  flow  boundary  conditions.  The 
remaining  set  of  equations  are  written  in  terms  of  (a,6,y)  coordinates  which 
feature  the  surface  located  along  an  a  =  constant  surface.  The  appropriate 
transformation  T3:  (a  ,3  ,y  )  +  (5  ,n  ,? )  as  a  function  of  edge  number  is: 


edge 

1: 

oc  =  5 

3  =  n 

Y  =  ? 

edge 

2: 

a  =  -n 

3  =  C 

Y  =  C 

edge 

3: 

a  =  -£ 

3  =  -n 

Y  =  4 

edge 

4: 

a  =  n 

3  =  -5 

Y  =  C 
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At  surface  points  the  quantities  are  P  =  £n(p),  V 
(entropy)  are  advanced  using  the  relations: 

3S  _  S  3S 
3y  w  36 


a  v  -  a  u  and  s 

y  x 


-  4  c  -  *!  -  «2  ♦  «j  -  ^ 

PS 


(12) 


where 


B  =  6  w  +  6  u  +  6  v 

z  a  y 


V  (az* 2  +  a*  +  ay2 *)w2  -  ( a*  +  ay2)a2 


Form  1 


A0  =  {"2az  +  a^^w2  "  ^ 


A1  =  +  a**)/(w2  "  ^ 

3a  3a  3a 

61  =  pw  [W  _  +  u  _  +  v  —  - 


*2=pWAl  ^az  fa  +  ax  fa  +  ay  3 a  ] 
*3  =  (“zA0  +  “x2  +  «/>£ 


3a 


3a 


64  =  -U  5F  +  V  9-T 


.  _  1  r  1  3p  Qi  3u  3V\  -1 

hl  "  pw  *-  J  36  pB^ay  36  ”  ax  36  ■* 


3  p  .  r  „  dU  ,  r  dV  dW  1 

h2  "  S0  36  P  ‘  S1  36  S2  36*  3  36  ^ 


3  u 


3  V 


3W 


s0 


-  BzA0  -  axBx  -  Vy 


S1  =  wBxAl  "  “xB 


s2  -  WBjAj  -  c,yB 


(13a) 


s3  =  wezAi  '  BAo 
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Form  2  ,  9 g,  9g?  9gd  3a 

hl  =^C  J  {V33  "ay  96'“+ax  83  }+p^y9F 


9a 

, 

’x  96 


}  ] 


(13b) 


h?  “  t  C  I-?  -  e  {  PW  (6/J)  +  PU  a|  (6  JO)  +  PV  ^  (3  V/J)  } 


where 


+  P  {  0’  gg  (S2/J),  gg  (3X/J),  gg  (By/J)}  -1 


? 

k,  „  w  k.  uwk 

t  =  J  {  (2  -  ~  q  )wAr  — -  A1  -  An , 


j  vwklAl 

O’  p  “  A1  ax’  "  ay  J 


e  {  l>w,u,v  }  ,  k^  ypj- 


p  =  ~Y-  ]T  (perfect  gas) 


As  indicated  above,  two  different  forms  are  available  for  evaluating  h^  and 
h2.  Both  forms  are  algebraically  equivalent;  however,  form  1  differences 
p,u,v,w  in  the  a  direction  while  form  2  differences  conservative  groupings  of 
these  quantities.  Form  1  is  found  to  be  more  robust  in  expansion  regions  and 
form  2  produces  superior  results  near  shocks. 


The  derivatives  which  appear  in  Equations  (12)  and  (13)  are  evaluated 
using  the  MacCormack  method.  In  the  a  direction,  first  order  accuracy  can  be 
obtained  using  one-sided  differencing  in  both  the  predictor  and  corrector 
step: 

predictor:  f  =  (f.x1  .  -  f.  .)/Aa 
a  1+1  ,J  i  >J 

corrector:  f  =  (f...,  .  -  f •  . )/Aa 
a  v  i+l, J  1,J 

Here  the  overbar  indicates  predicted  values.  It  second  order  accuracy  is 
required,  the  corrector  step  is  modified  as  follows: 


0% 


i+1  ,j 


-  f,  j)/4« 


(C 


1  ,J 


-  2f • 


i+1,  j 


(14) 


At  the  end  of  each  computational  step,  it  is  necessary  to  decode  the 
advanced  quantities  to  determine  p,p,u,v,w,  consistent  with  the  total  enthalpy 
constraint,  Equation  (2),  and  tangent  flow  boundary  condition.  The  required 
relations  are  given  by: 


p  =  exp(P) 

p  =  {exp[P-s(Y-l)]}1/Y 


w 


u 


•t/2(H0  -  Y  p/[(y-l  )p  ])  (ax2  +  ay2)-V2 

r  2~ 7  2-1 

La x  +  ay  +  az  J 
2  2 

(ayV  -  axazw)/(ax  +  ) 


(15) 
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?  ? 

V  =  -(a  V  +  a  a  w)/(a  +  a  ) 

'x  y  z  7  '  x  y7 

At  edge  points  lying  on  an  abutting  edge,  the  numerical  scheme  must  be 
modified  to  allow  access  to  the  information  in  the  adjacent  zone.  This  is 
accomplished  using  the  same  extended  computational  space  method  applied  at 
interior  points  lying  on  abutting  edges. 

c.  Juncture  Points 


At  a  juncture  point  two  surfaces  intersect.  Both  internal  and  external 
junctures  are  considered;  the  former  has  an  included  angle  of  less  than  180° 
while  the  latter  has  an  included  angle  greater  than  180°.  The  general 
strategy  for  treating  corner  points  is  to  define  a  mapping:  T^:  (a,8,y)  + 

(£,n,c),  the  inverse  of  which  transforms  the  corner  surface  to  a  =  0.  The 
characteristic  relations  developed  for  surface  points  can  then  be  applied  to 
the  corner.  In  order  for  the  transformation  to  be  non-singular  at  the  corner, 
it  is  necessary  to  envision  the  corner  as  slightly  rounded  on  a  sub-grid 
scale.  Since  the  transformation  needed  only  be  evaluated  at  discrete  grid 
points,  the  actual  corner  geometry  need  not  be  provided  and  only  the 
orientation  of  the  a  =  0  surface  at  the  corner  itself  need  be  specified.  Here 
it  is  required  that  the  direction  Va  bisect  the  corner  angle  in  the  crossflow 
pi ane. 


To  demonstrate  construction  of  T3,  consider*the  corner  1  coordinate 
system  which  is  illustrated  in  Figure  7.  In  the  neighborhood  of  the  corner 
point  a,  6,  y ,  are  described  by: 


a 


+  r  ?  +  (l-n)b 
"  L  (Anb  +  A?) 


]  ;  6 


=  JL  +  (n-l) . 

A?  An  ’ 


Y  =  A 


(16) 


Here  the  +,-  signs  apply  to  interior  and  exterior  corners  respectively.  The 
value  of  b  is  determined  by  requiring  that  Va  bisect  the  corner  angle  in 
(x,y,z)  space.  Noting  that  a  vector  bisecting  the  corner  angle  in  (x,y,z) 
space  is  given  by: 


VC  Vn 

I  vc  I  "  fvn 

and  that 


(17) 


i-Va  =  (a^Cx  +  arinx),  J*Va  =  (a^C  +  a^ny) 

where : 


=  l/(bAn  +  AC) 

a  =  -b/ ( bAn  +  AC) , 
n 


leads  to: 


?x  - 

bn 

X 

|  Vn  | 

-  n  I  vd 
x  1  ' 

S' 

bn 

y 

?y 

|  Vn  | 

0 

>Y 

1 

vc 

vn 


(18) 


Evaluating  at  corner  1  yields. 
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where:  m 


a 


a 


n 


=  m/(A5  +  bAn)  ,  g  =  -n/A£ 
=  nb/(A?  +  bAn )  ,  6^  =  m/An 
n  =  1 


(19) 


In  general  in,n  take  on  the  values  of: 


+1  corners  1  and  4 


+1  corners  3  and  4 


m 


n 


-1  corners  2  and  3 


-1  corners  1  and  2 


To  apply  the  MacCormack  scheme,  it  is  also  necessary  to  determine  T3  at  points 
adjacent  to  the  corner  (e.g.,  A,B  in  Figure  7).  Here  a  and  8  are  assumed 
parallel  to  £  or  n  and  T3  is  easily  determined.  For  example,  at  corner  1  the 
metrics  are: 


A:  a  =  ±  I/A£  ,  8  =  1/An ,  a  =  a  =  (3  =  8  =0 

C>  >1  '  I  S»  S 

3:  =  +  1/An  ,  =  1/AC  ,  =  8^  =  &  =  0, 


(20) 


where  the  top  sign  applies  to  interior  corners. 

In  applying  the  characteri Stic  relations  to  corner  points,  it  is  possible 
to  choose  from  forms  1  or  2,  and  first  or  second  order  differencing.  Also, 
corner  points  can  be  modeled  as  sharp  or  rounded.  At  a  sharp  corner,  the 
velocity  vector  is  aligned  with  the  corner  direction,  while  at  a  rounded 
corner,  arbitrary  flow  direction  is  permitted.  Application  of  the  surface 
characteri Stic  relations  simulates  a  rounded  corner.  Sharp  corner  modeling 
can  be  achieved  by  neglecting  the  calculated  V  and  substituting  the  constraint 
that  the  final  velocity  vector  be  parallel  with  the  corner  direction.  This 
implies  that  the  entropy  along  the  corner  is  constant.  Experience  suggests 
that  the  most  robust  corner  treatment  can  be  achieved  using  form  1,  first 
order  differencing,  and  the  rounded  corner  model. 

d .  Shock  Points 

The  shock  fitting  procedure  is  only  applicable  in  cylindrical 
coordinates.  It  is  designed  to  provide  an  outer  boundary  for  the 
computational  domain,  not  to  fit  imbedded  shocks.  Accordingly,  shock  points 
are  limited  to  edge  3  and  uniform  free  stream  properties  are  requi red  upstream 
of  the  shock.  At  the  shock,  flow  field  properties  as  well  as  the  shock 
location  and  slope  are  unknown.  The  correct  boundary  conditions  are  provided 
by  the  Ranki ne-Hugoni ot  conditions  which  relate  the  free-strearn  properties, 
the  shock  slopes,  and  the  properties  behind  the  shock.  An  analysis  of  the 
characteristics  associated  with  Equation  (1)  indicates  that  there  is  one 
admissable  characteristic  relation  at  the  shock.  This  relation,  when 
combined  with  the  Rankine-Hugoniot  relations,  results  in  a  system  of  equations 
to  advance  c  and  cz: 
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9cz  1  C2  r  Ccf,  3c  3c<fr  ,  .  D  C3 

‘  "  7  .  t  T  a7  “  »T"  J  +  K  “  r" ”  + 


3£  C  Cj  1  c  3?  3? 


Cj  V4 


(21) 


Evaluation  of  the  terms  appearing  on  the  riyht-hand-side  of  this  equation 
proceeds  as  follows: 

eY  =  (cos<j>  +  ^  si ntf> ) / v 

ey  =  ( si n<f»  -  ^  cos<j))/v 


ez  =  “cz,/v 


■  (  1  +  <c-)  +  czl 


}k 


(22) 


=  (1-y)  p  =  (1-y)p 
.2 


kl  '  Y  P 


U  =  exu  +  eyv  +  ezw 

U=eu  +  ev  +  ew  U 
co  x  °°  y  °°  z  °° 


^  ”V(U-we  )  +  (1-e  ^)(w^-a^) 


An  = 


4/(y  +  1) 

II  i|2  k  i 

(  -  -  1)  [1  +4  +  --  (P  -  P»)  1/(1 
U-  a^  p2 


P/Pco  <  1-02 

P/P„  >  1.02 


*0  =  p 


C  {  (e  -  3*)U  -  w}  An  - 


w  -  3*  (—  -  1)  U] 
U 


s0  =  WjT  +  (p-pj  +  PW< 
Si  =  u^ir  +  puw 

s2  =  v00ir  +  Pvw 
*  _  -yw  + 


0 


2  2 
w  -  a 


U  -  we. 


A1  =  wAq  +  y 


Z  (  ,  ^2 i  ^2 ’  ^4 1 
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*  M 

ll  ~  A1  p" 

★ 

wk,A,  * 

„  11  * 

2  p  o 

* 

llk^Aj 

£3  =  ~  ex 

* 

vk. 

£4  =  “T“  "  ey 


.  ap  ar  3 ( a  /J ) 

r  =  -j(  t-  [  <£  +  !!>  - « (-4—  * 


8a 


3(3  /J) 

— h—  > ] 


Co  ezsO  +  exsl  +  eys2 

C1  =  v  (_S0  +  ez  C0} 
r 

bo  1 

c~  =  -£=-v  [  C0(eycos<|>  -  ex  sin<(»)  +  s1  si -  s2  cos«j>] 


C3  _  exS2"  eySl 


_  j-  5z(sin<f,  ny  +  cos*  nx)  -  ny(sin<}>  5y  +  cos<j>  ?x)  -j 
- - (TJ1 


The  shock  calculation  procedure  provides  c,cz,  and  c  ,  while  the  user 
must  provide  <j>2,  <j>„  ,  <|>~  and  <t>2  .  To  calculate  metric  quantities,  the 
derivatives  c7,  c.  ,^c,.,  c77  anci  c,_  are  necessary.  These  are  determined  as 
follows:  Z  +  ++  ZZ  *z 


8C  1 

(j)  8n  <f> 


(cj+l  ~  cj-l*  1 
2An  <j> 


where. 


1 

8cz  _  « 

ZJ+1  '  ' 

CZ(f) 

*n 

8n 

2  An 

r  = 

1 

9  ca 

- (A) 

8n  \ 

1 

W 

♦ 

n 

2An<t>, 

c  = 

c 

-  c  , , a0 

-  C .  do 

ZZ 

zc 

M  2 

(j>z  3 

a2  = 

z’  a3  =  ^ 

z 

z . 


nj 


nj  n 


(23) 
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<l>n  =  ( 4>2  “  ^3 ) G ‘  (n) 


At  zone  boundaries,  the  mesh  may  be  discontinuous  and  the  above 
expressions  could  yield  different  c  derivative  values  at  abuttiny  shock 
points.  At  such  locations,  the  procedure  for  the  calculating  c  and  c  is 
modified.  To  illustrate  the  method,  consider  the  following  four  succe^iive 
shock  points;  j'-l,j\  j",  j"-l  where  j'  and  j"  abut.  The  '  points  are 
located  in  zone  I  and  “  points  in  zone  II.  c.  and  c  are  calculated  as 
follows:  ♦ 


c  =  c,  +  (c,  -  c,  )f 


4  -i  1  4 


•J"  TJ‘  yR  VL  VR 


4i  4d  c 


(24) 


C 4 z  - ,  C4ZD  +  "  c<j>_  )fc 


zj 


ZL  ZR 


where: 


(c..|  C  ♦  1  ■*  )  /  Ar|  <j>  I  >  C  ( ^  i  *'4.1  C  -  H  )  /  Ar)  <j>  1 1 

<J>L  J  J  n j  *  9r  J  j  n j  >■ 


C*  =  (  C  7  -  C7  )/AT).  I  ’  CA  =  (C7  "  C7 

9-,  z 1  z-,«  i  9  ^  9-,  Z  h,,  Z-„ 

ZL  J  J  _1  dj ,  rzR  j  +1  j 

f'c  =  An 1 1 4^  1 1  j  /  ( An  ^  4>  j  +  An^4  n  ) 

n  J‘  n  j" 


-  c7  )/AnU4  II 

V 


On  the  leeside  of  a  body,  the  shock  becomes  a  Mach  line  and  the  shock 
tracking  procedure  is  replaced  by  a  Mach  line  trackiny  algorithm.  Mach  line 
tracking  is  accomplished  by  requiring  that  the  tree-stream  Mach  Number  normal 
to  the  shock  be  unity,  or: 


q 


00 


•  n  = 


(25) 


where:  n=e  -  ^  e  -ce 

r  c  <j>  z  z 

Using  known  values  of  c,  c^  and  solving  for  cz  produces: 


K-  E  VJ  +  -c-  vj2  +  (w 2  -  a_2)(l  +  cj/c2) 


[w  2  -  a  2) 


(26) 


e.  Symmetry  Plane 

A  symmetry  plane  is  any  plane  about  which  the  dependent  variables  are  odd 
or  even  functions.  This  allows  properties  to  be  constructed  on  a  fringe 
plane,  located  outside  of  the  computational  domain,  and  thus  permits  symmetry 
plane  points  to  be  advanced  using  the  standard  MacCormack  algorithm.  To 
simplify  the  computational  scheme,  the  symmetry  plane  is  required  to  satisfy 
the  following  constraints: 
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1.  The  symmetry  plane  normal,  n,  must  have  a  null  z  component. 

2.  On  the  symmetry  plane  Va*V3  =  0. 

3.  Either  Va  or  V3  must  be  parallel  to  n. 

These  restrictions  allow  metric  values  to  be  determined  on  the  fringe  plane 
using  odd  or  even  constructions. 

To  construct  fringe  plajie  values  of  the  dependent  variables,  two  cases 
must  be  considered:  (1)  Va*n  =  0,  (2)  V3*n  =  0.  Here,  F  is  a  vector  normal 
to  the  symmetry  plane.  In  the  first  case,  the  symmetry  plane  is  parallel  to 
a  3  =  constant  plane  while  in  the  second,  it  is  parallel  to  an  a  =  constant 
plane.  To  illustrate  the  procedure  for  constructing  fringe  plane  properties, 
let  point  1  be  on  the  symmetry  plane,  2  be  the  adjacent  point  inside  the 
computational  domain,  and  0  be  the  fringe  plane  point.  In  either  case,  p  and 
p  are  even  functions  leading  to: 


P0  ?2 

The  velocity,  q,  follows  by  noting  that  qp  and  q.  are  odd  and  even  functions, 
where  the  subscript  n  or  t  refers  to  the  vector  component  normal  and  tangent 
to  the  symmetry  plane.  If  the  symmetry  plane  satisfies  Vcwt  =  0,  then  (va) 
and  (V3)n  are  odd  while  (va)n  and  (v3)t  are  even*.  For  V3*n  =  0,  the  odd  ana 
even  behavior  of  the  metric  normal  andLtangenti al  components  is  reversed. 
These  considerations  lead  to  the  following  scheme  for  constructing  fringe 
plane  velocity  and  metric  values: 

Case  (1):  3  =  constant  parallel  to  the  symmetry  plane,  (i.e.,  Va*F  =  0) 


.  2(Vcu*V31) 

Vctn  =  Va?  -  [  - ?  ■  ]V3, 

v$i 

2(V32*  VM 

V3q  =  [  - 7> -  ]V3,  -  V3o 

V3j 

t  _  t  r  o  (V3r^2} 

%  Qo~[2  /p  ]v$-i 

76x 

Case  (2):  a  =  constant  parallel  to  the  symmetry  plane. 
2(Vct2*Va1) 

VcIq  =  [ - 2 -  ]  4a^  -  Va^ 

Val 

2  ( V3n  * Va. ) 

V3q  =  V32  -  [  -  -2 — —  ]  Vox 

V«1 


(28) 


(i.e.,  V3*n  =  0) 


(29) 
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r  2(Varq0}  1 

q0  =  q2  -  [ - 1 - JVa 

Va^ 


1 


The  MacCormack  scheme  is  not  symmetric  and  direct  application  of  this 
algorithm  to  symmetry  points  can  lead  to  an  unstable  calculation.  On  the 
symmetry  plane,  the  algorithm  is  altered  at  both  the  end  of  the  predictor  and 
correct  step  to  ensure  that  the  decoded  velocity  component  normal  to  the 
symmetry  plane  is  zero.  This  is  accomplished  by  redefining  the  third  and 
fourth  components  of  the  U  vector,  U3  and  u^ ,  as  follows: 

u3'  =  (Va  ’  ayV/(Vy  "  ayex} 


=  (a  V  -  B  V  )/(a  B  -  a  8  ) 
1  x  e  xa/vxy  yx' 


where: 


r  lUa  + 

=  {  3  x 

u4ay  if 

symmetry 

pi  ane 

is 

paral lei 

to 

Va 

0 

if 

symmetry 

pi  ane 

i  s 

paral 1  el 

to 

VB 

- 1  0 

if 

symmetry 

pi  ane 

is 

paral 1  el 

to 

VB 

1  U3Bx  + 

U .  3  if 
4  y 

symmetry 

pi  ane 

is 

paral 1  el 

to 

Va 

(30) 


here  (U31  ,  U41)  and  (U3,  U4)  are  new  and  old  values  respectively. 


When  treating  axisymmetric  problems,  it  is  desirable  that  in  cylindrical 
coordinates,  all  of  the  dependent  variables  be  identical  on  each  constant  <j> 
surface.  To  accomplish  this,  a  different  procedure  for  determining  fringe 
plane  circumferential  and  radial  velocity  components  is  needed: 


q0  (q2  "  qV 


2 


[VBj- (q2  -  q1)]VB1  + 


VB 


2 

1 


ql 


(31) 


Application  of  the  above  removes  the  necessity  of  setting  the  symmetry  plane 
velocities  to  zero. 


4.  SPECIAL  PROCEDURES 

In  order  to  complete  calculations  on  complicated  configurations,  it  is 
necessary  to  introduce  a  number  of  special  procedures.  These  deal  with  the 
treatment  of  surface  edges  and  application  of  artificial  viscosity.  Surface 
edges  can  occur  at  a  cowl  lip,  wing  edge,  or  as  a  result  of  a  surface  slope 
discontinuity.  The  properties  downstream  of  an  edge  are  estimated  using  a 
local  analysis,  and  special  differencing  is  applied  to  promote  robustness. 

a.  Surface  Edge  Analysis 


As  a  calculation  is  marched  down  the  axis  of  a  missile,  new  wing  or  cowl 
surfaces  may  be  encountered.  The  computational  zones  must  be  chosen  to  ensure 
that  these  surfaces  lie  along  zone  boundaries.  Along  edge  boundaries,  two 
points  are  associated  with  each  point  in  computational  space,  one  with  each 
zone.  These  two  sets  of  points  are  used  to  describe  the  two  sides  of 
surfaces.  New  surfaces  can  be  encountered  anywhere  during  a  calculation,  and 
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points  which  at  one  step  are  interior  points,  may  in  the  next  step  become 
surface  points  and  vice-versa.  A  leading  edge  point  is  a  surface  point  which 
in  the  previous  step,  was  an  interior  point  while  a  trailing  edge  point  is  an 
interior  point  which  in  the  preceding  step  was  a  surface  point.  In  addition, 
surface  slope  discontinuities  can  be  encountered  during  the  course  of  a 
cal cul ation. 

A  local  analysis  is  generally  applied  at  surface  edges  which  models  the 
flow  in  the  vincinity  of  the  edge.  When  the  flow  normal  to  the  wing  edge  is 
sufficiently  supersonic,  this  analysis  is  based  on  the  oblique  shock  or 
Prandtly-Meyer  expansion  relations.  In  other  cases,  an  empirical  treatment  is 
used.  Application  of  the  local  analysis  improves  robustness  and  accuracy  near 
surface  edges.  The  computational  algorithm  proceeds  by  completing  the  step  in 
which  the  surface  edge  is  encountered  without  taking  the  edge  into  account. 
During  this  step,  old  values  of  the  edge  derivatives  are  used  in  conjunction 
with  the  updated  edge  locations.  The  computed  flowfield  properties  at  the  end 
of  the  step  are  taken  as  the  conditions  immediately  upstream  of  the  surface 
edge  and  the  local  analysis  is  applied  to  turn  the  velocity  vector  parallel  to 
the  downstream  surface  slope. 

The  first  step  of  the  edge  analysis  is  to  calculate  vectors  normal  to  the 
surface  upstream  jand  downstream  of  the  edge  (i.e.  n+,n-  in  Figure  8).  At  the 
leading  edge,  n+,7F-  are  constructed  from: 

~n-  =  +  (Tx?i-)  (32) 

n+  =  (Va)+ 


The  sign  associated  with  n-  is  chosen  to  insure  n+,n-  >  0,  where  x  is  a  vector 
tangent  to  the  edg_e.  Here  (Va)  is  evaluated  downstream  of  the  edge.  At  the 
trailing  edge,  n+,n-  are  evaluated  using  the  values  of  Va  upstream  and 
downstream  of  the  surface  edge  respectively.  n+  is  perpendicular  to  the 
downstream  zone  edge  which  approximately  bisects  the  trailing  edge  angle.  At 
surface  discontinuities,  n+,n-  are  evaluated  using  upstream  and  downstream 
val ues  of  Va. 

After  determining  n+  and  n- ,  the  turn  angle,  o,  and  the  velocity 
components  normal  (qn)  and  tangent  (qt)  to  the  edge  are  calculated  from: 


o  =  cos"1  [ln-+--Q:->] 
I n+  I  I  n 


qn_  =  (s_  •  q)/  |s_| 

(t  •  q) / | x | 

where .  "x  =  (  n+)x(  n-) 


s  =  xxn- 


(33) 
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Using  the  known  upstream  values  of  p,  qn_  and  p,  the  local  analysis  procedure 
of  Refs.  9,  11  provides  downstream  values  of  these  quantities.  The  final 
downstream  velocity  vector  is  constructed  from: 


where: 


s+  =  xxn+ 


(34) 


At  the  trailing  edge,  the  above  procedure  is  applied  only  if  the  flow 
normal  to  the  trailing  edge  on  both  surfaces  is  supersonic.  Otherwise, 
properties  downstream  from  the  trailing  edge  are  determined  by  averaging 
adjacent  interior  point  properties.  When  both  surfaces  associated  with  a 
trailing  edge  feature  supersonic  flow  normal  to  the  edge,  the  shock  or 
expansion  turning  procedure  is  applied  to  each  surface  and  results  are 
averaged  to  determine  final  properties  downstream  of  the  trailing  edge. 


A  different  procedure  is  required  to  turn  the  flow  at  a  corner  leading 
edge  point.  In  the  case  of  an  external  corner,  a  method  similar  to  that 
applied  at  a  regular  surface  edge  point  is  used,  but  the  normal  vector  is 
taken  in  a  direction  bisecting  the  corner  angle.  At  an  internal  corner,  a  two 
step  procedure  is  used  first  to  turn  the  flow  into  a  plane  perpendicular 
to  V£  ,  and  then  to  rotate  it  in  this  plane  until  it  is  also  normal  to  Vn.  The 
required  turn  angle  for  aligning  the  velocity  perpendicular  to  vc  is  given  by: 


0 


±  sin 


(vc-q) 

M|q| 


(35) 


Here  the  plus  sign  applies  at  corners  1  ,< 
elsewhere.  This  initial  turn  produces  a 

s  =  q  -  larlli  vc 

|VC  | 

The  necessary  turn  required  to  align  the 
direction  is  given  by  : 

0  =  cos“1[s1*  l2/(  ITJ  )T2 1 ) ] 

Here  S£  =  -(Vnxvc)  and  is  the  final  flow 
corners  1  and  2,  Vn  •  s,>  0  corresponds 
implies  an  expansion.  The  sign  conventi 


while  the  minus  sign  is  used 
final  flow  direction  of: 

(36) 


velocity  vectory  with  the  corner 


(37) 

direction.  At  junctions  located  on 
o  a  compression  and  Vn  •  S,  <U 
>n  is  reversed  on  corners  3and  4. 


At  a  corner  trailing  edge  point  either  one  or  two  surfaces  may 
disappear.  If  one  surface  disappears,  the  velocity  vectors  on  each  side  of 
the  vanishing  surface  are  alligned  with  the  downstream  zone  corner  direction, 
given  by  Vnxvc.  This  is  accomplished  using  the  two  step  leading  edge  turn 
procedure  described  above.  The  trailing  edge  downstream  conditions  are 
obtained  by  averaging  the  flow  properties  occuring  on  each  side  of  the 
vanishing  surface.  When  two  surfaces  disappear  simultaneously,  the  standard 
trailing  edge  point  procedure  described  above  is  applied.  However,  the 
surface  normal  vector  is  taken  in  the  direction  bisecting  the  corner  angle. 
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b.  Surface  Edge  Differencing 

The  presence  of  surface  edges  requires  the  introduction  of  special 
differencing  procedures  for  advancing  points  labeled  T  and  I  in  Figure  9.  In 
all  cases,  modified  procedures  must  be  applied  at  points  labeled  I,  which  have 
two  neighboring  T  points  on  the  upper  and  lower  surfaces  of  the  wing.  I 
points  are  advanced  twice  using  each  of  the  neighboring  points,  and  following 
the  corrector  step  the  two  resulting  sets  of  property  values  are  averaged.  In 
situations  where  strong  shocks  or  expansions  are  attached  to  surface  edges, 
differencing  is  modified  in  the  following  fashion: 

(1)  Alteration  of  differences  along  zone  boundaries.  Both  T  and  I  points 

are  of  concern  here.  At  the  surface  edge,  a  strong  shock  or  Prandtl-Meyer 
expansion  may  occur.  Under  these  conditions,  T  points  are  not  effected  by  I 
points  and  vice-versa.  This  correct  domain  of  dependence  can  be  imposed  at 
point  I  by  applying  a  one-sided  difference  along  the  zone  boundary  directed 
away  from  the  surface  edge.  At  point  T,  this  procedure  leads  to  an 
instability  and  differences  calculated  using  points  T  and  I  are  set  to  zero. 

(2)  Normal  surface  derivative  damping  at  leading  edges.  Large  property 

jumps  can  occur  at  the  leading  edge  due  to  the  existence  of  a  shock  or 

expansion.  The  surface  normal  derivatives  at  T  type  points  are 

unrealistically  large,  reflecting  the  fact  that  the  difference  is  being  taken 
across  a  discontinuity.  Use  of  calculated  normal  derivative  values  in  the 
vicinity  of  the  leading  edge  leads  to  an  inaccurate  solution.  More  realistic 
values  can  be  obtained  by  damping  the  normal  derivatives  near  the  leading  edge 
(i.e.,  reducing  the  derivative  by  an  empirically  determined  corrective 
factor).  The  effect  of  normal  derivative  damping  is  illustrated  in  Reference 
12.  The  user  is  required  to  select  the  number  of  step,  K,  for  which  the 
derivatives  will  be  altered.  The  damping  factor  is  calculated  from: 


F k  =  [  MAX{1  -  0}]2  (38) 

where  k  is  the  number  of  step  following  the  occurrence  of  the  jump, 
c.  Artificial  Viscosity 


In  computations  featuring  separation  or  highly  swept  wings,  large  vortex 
structures  developed  in  the  flow  field.  In  such  circumstances,  it  is  often 
necessary  to  introduce  an  artificial  viscosity  into  the  calculation.  This 
smooths  the  gradients  surrounding  such  phenomena  and  leads  to  a  more  robust 
computational  procedure.  Smoothing  can  be  applied  to  all  points  except  shock 
points.  In  each  case,  the  advanced  quantities  are  smoothed  following  the 
corrector  step.  At  interior  point,  a  Schuman  filter  with  a  density  switch  is 
implemented  as  follows: 


=  tic 


+  [0C  . 
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(39) 
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where: 


is  the  advanced  flux 
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(40) 


and  Cv,  Cu  are  smoothing  factors  specified  by  the  user. 

*  J 

Wall  points  are  smoothed  using  a  weighted  average  of  neighboring  wal 


poi nts : 

t.  = 

J 

A0C  + 
J 

+ 

cOc  . 

J-1 

where: 

A  =  1 

-B-C 

B  =  | 

S+i 

(°j+l  - 

pj)/(pJ+l 

+  pj)  | 

c  -  1 

S-i 

(pj  ■  °0 

-l)/(pj  + 

PM>I 

In  the 

j-l»j+l 

above, 
poi nts 

CWJ  is  a  user 
are  depicted 

specified 
in  Figure 

constant 
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(41) 


(48) 


the  neighboring  points  are  not  both  wall  points  and  the  above  procedure  is 
altered.  Figure  10c- and  lOd  indicate  the  points  used  in  applying  smoothing 
which  is  implimented  by: 


tL 


5.  RESULTS 


C®j?  +  Cw.<^+^  +  ^)]/(1  +  3Cw.> 
J  J 


(43) 


This  chapter  presents  cases  which  have  been  run  with  the  multi  zone  code 
(MUSE).  All  computations  were  performed  on  a  cyber  170/720  and  required 
approximately  5(10-3)  seconds  per  point  per  step.  The  initial  data  plane  is 
generated  using  the  approximate  conical  solver  of  Reference  16.  The  actual 
body  geometry  used  as  well  as  the  input  data  list  are  provided  in  Reference 
13.  The  mesh  size  employed  for  each  zone  is  indicated  by  (NXM)  where  N  and  M 
are  the  number  of  points  in  the  radial  and  circumferential  directions 
respectively.  Unless  otherwise  indicated,  smoothing  was  not  applied,  and  one¬ 
sided  differencing  was  used  at  all  fin  tips. 

The  normal  force  and  center  of  pressure  computations  for  an  airplane  type 
configuration  at  Mach  2,  are  compared  to  the  experimental  data  of  Reference  17 
in  Figure  11.  This  case  features  attached  shocks  on  the  wing  and  tail 
surfaces.  The  calculation  was  carried  out  using  a  single  (15x15)  zone 
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forwards  of  the  wing  and  two  (29x20)  zones  over  the  remainder  of  the 
configuration.  To  prevent  the  flow  from  becoming  subsonic  at  the  body-wing 
junction,  the  leading  edge  turning  angle  at  this  location  was  damped.  The 
computed  and  measured  normal  force  coefficient  and  center  of  pressure  for  the 
body-wing  and  body-wing-tail  configurations  at  an  incidences  of  5°,  10°  and 
15°  are  in  good  agreement.  The  computed  crossflow  velocities  and  pressures  at 
the  wing  and  tail  trailing  edges  at  an  incidence  of  10°  are  shown  in  Figure 
12. 


In  Reference  18,  a  tangent  ogive  body  with  cruciform  delta  fins  was 
tested  at  Mach  3.7.  This  configuration,  which  is  shown  in  Figure  13,  features 
attached  shocks  at  the  wing  leading  edges.  Numerical  results  are  compared  to 
experiment  in  Figure  13  for  the  case  of  the  model  in  a  plus  roll  orientation 
and  an  incidence  of  7.8°.  A  single  (15x15)  zone  was  applied  forward  of  the 
wing  and  two  (40x19)  zones  were  employed  over  the  remainder  of  the  body.  All 
fin  surfaces  show  reasonable  agreement  between  calculation  and  experiment, 
although  problem  areas  occur  near  the  body-tail  juncture.  This  area  is 
strongly  influenced  by  viscous  effects  and  lack  of  close  agreement  is  not 
surprising.  The  computed  flow  field  at  an  axial  station  near  the  center  of 
the  wing  is  shown  in  Figure  14.  Clearly  evident  is  the  interaction  between 
the  shocks  generated  by  the  horizontal  and  windward  fin. 

Reference  19  provides  body  and  wing  pressure  data  on  a  cruciform  delta 
wing  configuration  at  Mach  2.7.  Calculations  were  performed  on  this  geometry, 
which  is  illustrated  in  Figure  15,  at  an  incidence  of  10°,  with  the  model  in 
the  plus  roll  orientation.  Under  these  conditions,  shocks  at  the  wing  leading 
edges  were  slightly  detached.  A  uniform  (15x15)  zone  was  used  forward  of  the 
wing  and  two  (40x19)  zones  were  applied  for  the  remainder  of  the 
calculation.  Comparison  of  fin  and  body  surface  pressures  are  shown  in 
Figures  15  and  16  respectively.  Both  body  and  wing  surface  pressures  are  in 
reasonable  agreement  with  experiment.  However,  the  body  surface  pressure 
results  in  under-predicting  the  rate  at  which  fin  disturbances  propagates  over 
the  body.  The  crossflow  velocity  and  pressures  near  the  wing-center  are 
exhibited  in  Figure  17. 

Calculations  have  been  performed  for  the  two  swept  wing  configurations, 
tested  in  Reference  20  and  shown  in  Figure  18.  These  models  were  run  at  an 
incidence  of  6°,  Mach  numbers  of  2.5  and  4.5,  and  feature  the  same  planform 
but  different  fin  thickness.  Such  conditions  result  in  detached  leading  edye 
shocks  in  all  cases.  A  (25x25)  mesh  was  applied  forwards  of  the  wing  and  two 
(29x20)  zones  were  used  for  the  remainder  of  the  body.  Calculated  and 
measured  wing  surface  pressures  are  shown  in  Figure  18  and  agree  reasonably 
well  with  experiment.  However,  near  the  wing  leading  edge,  computed  values 
are  larger  than  measured  ones.  Figure  19  provides  measured  and  calculated 
surface  pressure  on  the  windward  and  leeward  side  of  the  body.  Upstream  of 
the  wing,  calculation  and  experiment  are  in  good  agreement,  but  the  calculated 
jump  in  body  pressure  on  the  windward  side  caused  by  the  presence  of  the  wing 
occurs  further  aft  than  the  measured  one,  particularly  at  Mach  4.5. 

The  crossflow  velocity  and  pressure  field  is  illustrated  in  Figures  20 
and  21  at  several  axial  stations  along  the  wing  section  for  thick  fin  models 
at  Mach  numbers  of  4.5  and  2.5,  respectively.  At  the  axial  station,  z  =  56, 
in  each  figure,  the  formation  of  a  detached  wing  shock  is  visible  below  the 
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wing.  Axial  locations  z  =  62  and  z  =  65  illustrate  the  formation  of  an 
expansion  associated  with  the  wing  surface  slope  discontinuity.  A  decaying 
shock  wave  is  visible  in  the  flow  field  at  a  station  near  the  wing  trailing 
edge,  at  z  =  75.  In  the  higher  Mach  number  case  the  shock  is  positioned 
nearer  to  the  wing  surface.  This  produces  the  strong  wing  surface  pressure 
gradients  which  are  visible  at  this  Mach  number  in  Figure  18. 

Figure  22  illustrates  a  swept  wing  configuration  with  a  vertical  tail 
located  on  the  wing.  Calculations  were  carried  out  at  Mach  2.86  at  incidences 
of  4.05°  and  8.47°.  Under  these  conditions,  the  wing  leading  edge  shock  was 
detached.  This  calculation  was  carried  out  using  a  single  (18x19)  zone  up  to 
the  wing.  Two  (40x19)  zones  were  applied  over  the  forward  section  of  the  wing 
and  four  zones,  two  (30x19)  and  two  (15x19)  zones  were  employed  over  the 
vertical  tail  section.  The  tail  thickness  was  neglected  in  order  to  simplify 
the  geometry  description.  The  calculated  normal  force  coefficient  and  center 
of  the  pressure  are  shown  in  Figure  22  and  agree  well  with  experiment. 
Crossflow  velocities  and  pressures  at  an  axial  station  slightly  forward  of  the 
tail  are  shown  in  Figure  23  along  with  those  for  a  similar  configuration 
without  vertical  tail.  The  presence  of  the  tail  is  seen  to  produce  large 
changes  in  the  leeward  flow  field.  A  detached  shock  can  be  seen  lying  below 
the  outboard  wing  section  and  a  crossflow  shock  occurs  on  the  leeside  of  the 
wing  in  both  cases. 

Results  for  a  missile  with  wrap-around  fins  at  Mach  2.5  and  3  are 
presented  in  Figure  24  to  26.  The  missile  geometry  is  shown  in  Figure  24 
along  with  the  measured  normal  force  coefficient  on  the  spinning  model  of 
Reference  22.  A  single  (15x25)  zone  was  employed  forward  of  the  wing  while 
four  (40x10)  zones  were  applied  over  the  fins.  Mesh  clustering  was  used  in 
the  radial  direction  to  increase  the  number  of  points  in  the  vicinity  of  the 
fins.  As  is  indicated  in  Figure  24,  good  agreement  is  obtained  between  the 
computed  and  measured  normal  force.  A  discrepancy  between  calculated  and 
measured  center  of  pressure  is  possibly  attributable  to  the  fact  that  the 
experimental  model  was  spinning.  The  rolling  moment,  which  was  measured  in 
Reference  23  is  exhibited  in  Figure  25  and  does  not  agree  quantitatively  with 
calculation.  Flowever,  results  correctly  predict  trends  with  Mach  number  and 
incidence.  Crossflow  velocities  and  pressures  in  the  vicinity  of  the  wrap¬ 
around  fin  trailing  edge  are  shown  in  Figure  26  at  a  Mach  number  of  3  and  an 
incidence  of  8  degrees. 

6.  CONCLUDING  COMMENTS 

The  MUSE  code  described  in  this  report  is  applicable  to  missile  shapes  in 
supersonic  flight.  The  explicit  MacCormack  scheme  is  used  to  advance  interior 
points  while  a  characteristic  boundary  condition  treatment  is  applied  at 
surfaces  and  the  shock.  A  local  analysis  is  invoked  at  wing  edges  and  at 
surface  slope  discontinuities,  which  leads  to  accurate  surface  pressure  values 
near  these  locations  and  adds  to  the  robustness  of  the  procedure.  Meshes  are 
generated  algebraically  and  complex  shapes  are  handled  by  dividing  the 
computational  domain  into  a  number  of  simple  zones.  MUSE  differs  from  SWINT12 
by  allowing  the  zone  shape  to  be  completely  specified  and  by  permitting  zones 
to  be  connected  in  a  relatively  arbitrary  manner.  The  MUSE  code  removes  the 
thin  fin  assumption  present  in  SWINT  and  is  capable  of  treating  a  broad  range 
of  configurations,  many  of  which  are  beyond  the  scope  of  SWINT.  Examples  of 
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such  shapes  are  illustrated  in  Figures  22  through  26  and  include  a  wrap-around 
fin  missile,  and  airplane  type  configurations  with  vertical  tails  on  the  wing. 
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FIGURE  1.  CARTESIAN  AND  CYLINDRICAL  COORDINATE  SYSTEM 
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FIGURE  2.  WRAP-AROUND  MESH 
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FIGURE  3.  MULTIPLE  ZONE  MESH 
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FIGURE  4.  GENERALIZED  QUADRILATERAL  ZONE  STRUCTURE 


27 


NSWC  TR  85-484 


CONTINUOUS 


DISCONTINUOUS 


PRESSURE  CONTOUR 


MESH 


MESH 


-  CONTINUOUS  MESH 

-  DISCONTINUOUS  MESH 


FIGURE  5.  COMPUTED  PRESSURE  DISTRIBUTION  ON  A  CONE  USING  BOTH  A  CONTINOUS 
AND  DISCONTINUOUS  MESH  (MACH  =  3,  CONE  ANGLE  =  7°,a=8°) 
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FIGURE  6.  CALCULATED  PRESSURE  ON  A  BICONIC  WITH  A  UNIFORM  AND  A  DISCONTINUOUS 
MESH  (MACH  =3,  CONE  ANGLE  =7.5°,  20\a=0°  ) 
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FIGURE  7.  JUNCTURE  COORDINATE  SYSTEM 
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FIGURE  9.  SPECIAL  POIMTS  ADJACENT  TO  EDGE  SURFACES 
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FIGURE  10.  POINTS  USED  IN  SURFACE  SMOOTHING 
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FIGURE  11.  CALCULATED  AND  MEASURED  NORMAL  FORCE  COEFFICIENT,  CN,  AND  CENTER 

OF  PRESSURE,  ZAC/L,  ON  THE  BODY-WING  AND  BODY-WING-TAI L  CONFIGURATION 
OF  REFERENCE  17,  AT  MACH  =  2  ANDa=  5°,  10°  AND  15° 
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FIGURE  12.  CALCULATED  CROSSFLOW  PLANE  VELOCITY  AND  PRESSURE  CONTOURS  AT  BOTH 
THE  TRAILING  EDGE  OF  THE  WING  AND  OF  THE  TAIL  FOR  THE  CONFIGURATIONS 
OF  FIGURE  11,  MACH  =  2,  G=  10° 
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FIGURE  13.  CALCULATED  AND  MEASURED  FIN  SURFACE  PRESSURES  ON  A  CLIPPED  DELTA 
CONFIGURATION, &  =  7.8°  #  MACH  -3.7.  DATA  IS  FROM  REFERENCE  18. 

CURVES  HAVE  A  ZERO  REFERENCE  SHIFTED  BY  ONE  FOR  EACH  SUCCESSIVE 
SPANWISE  LOCATION  .  THE  SYMBOLS  ARE  EXPERIMENTAL  DATA  AND  THE 
LINES  ARE  CALCULATION. 
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FIGURE  14.  CALCULATED  CROSSFLOW  VELOCITIES  AND  PRESSURE  CONTOURS  FOR  THE 
CLIPPED  DELTA  FIN  CONFIGURATION  OF  REFERENCE  18  IN  A  PLUS  ROLL 
ORIENTATION,  MACH  =  3.7,  AND  a  =  7.8° 
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FIGURE  15.  CALCULATED  AND  MEASURED  SURFACE  PRESSURES  ON  THE  CRUCIFORM  DELTA 
CONFIGURATION  OF  REFERENCE  19  IN  THE  PLUS  ROLL  ORIENTATION,  MACH  = 
2.7,  a  =  10°.  SUCCESSIVE  CURVES  HAVE  A  ZERO  SHIFT  OF  UNITY*  THE 
SYMBOLS  ARE  EXPERIMENTAL  DATA  AND  THE  LINES  ARE  CALCULATION. 
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FIGURE  16.  CALCULATED  AND  MEASURED  BODY  SURFACE  PRESSURES  ON  THE  CRUCIFORM 
DELTA  CONFIGURATION  OF  REFERENCE  19  IN  THE  PLUS  ROLL  ORIENTATION, 
MACH  =  2.7,  a=  10”.  THE  SYMBOLS  ARE  EXPERIMENTAL  DATA  AND  THE  LINES 
ARE  CALCULATION. 
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FIGURE  17.  CALCULATED  CROSSFLOW  VELOCITIES  AND  PRESSURE  CONTOURS  FOR  THE 

DELTA  FIN  CONFIGURATION  OF  REFERENCE  19  IN  A  PLUS  ROLL  ORIENTATION, 
AT  MACH  =  2.7,  a  =  10° 
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FIGURE  18.  CALCULATED  AND  MEASURED  FIN  SURFACE  PRESSURES  ON  THE  THIN  AND 
THICK  WING  CONFIGURATIONS  OF  REFERENCE  20,  AT  MACH  =  2.5,  4.5 
AND  a  =  6° 
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FIGURE  19.  CALCULATED  AND  MEASURED  BODY  SURFACE  PRESSURES  ON  THE  THIN  AND 
THICK  WING  CONFIGURATIONS  OF  REFERENCE  20,  AT  MACH  =  2.5,  2  4.5 
AND  a  =  6°.  THE  SYMBOLS  ARE  EXPERIMENTAL  DATA  AND  THE  LINES  AE  ARE 
CALCULATION. 
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FIGURE  20.  CALCULATED  CROSSFLOW  VELOCITIES  AMD  PRESSURE  CONTOURS  ON  THE  THICK 
WING  CONFIGURATION  OF  REFERENCE  20  AT  MACH  =  4.5, Cl  =  6° 
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FIGURE  21.  CALCULATED  CROSSFLOW  VELOCITIES  AND  PRESSURE  CONTOURS  ON  THE  THICK 
WING  CONFIGURATION  OF  REFERENCE  20  AT  MACH  =  2.5,a=  6° 
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FIGURE  22.  NORMAL  FORCE  COEFFICIENT,  CN,  AND  CENTER  OF  PRESSURE,  ZAC/L,  ON 
THE  BODY-WING-TAIL  CONFIGURATION  OF  REFERENCE  21  AT  MACH  =  2.86 
AT  a  =  4.05°,  8.47° .  THE  LINES  ARE  THE  EXPERIMENTAL  DATA  AND  THE 
SYMBOLS  ARE  CALCULATION. 
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FIGURE  23.  CROSSFLOW  VELOCITIES  AND  PRESSURE  CONTOURS  ON  THE  BODY  WING  AND 
BODY-WING-TAIL  CONFIGURATION  OF  REFERENCE  21  AT  MACH  =  2.86  AT 
a  =  8.47°  AT  AN  AXIAL  STATION  SLIGHTLY  FORWARD  OF  THE  WING 
TRAILING  EDGE 
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FIGURE  24. 
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CALCULATED  AND  MEASURED  NORMAL  FORCE  COEFFICIENT  ON  A  WRAP-AROUND 
FIN  CONFIGURATION  TESTED  IN  REFERENCE  22,  AT  MACH  2.5  AND  3. 

THE  MACH  2.5  CURVE  IS  OFFSETBY  .5.  THE  LINES  ARE  EXPERIMENTAL  DATA 
AND  THE  SYMBOLS  ARE  CALCULATION. 
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FIGURE  25.  CALCULATED  AND  MEASURED  ROLL  MOMENT  COEFFICIENT  ON  THE  WRAP-AROUND 
FIN  CONFIGURATION  TESTED  IN  REFERENCE  24.  MEASUREMENTS  ARE  TAKEN 
FROM  REFERENCE  23.  THE  SYMBOLS  ARE  EXPERIMENTAL  DATA  AND  THE  LINES 
ARE  CALCULATION. 
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APPENDIX  A 

ANALYTIC  DERIVATION  OF  METRICS 
The  metric  quantities  £ 


nx  ny<  nz,  appearing  in  Equation  (1) 


x  ’y  ’  z 

must  be  evaluated  at  every  point  in  the  calculation.  In  addition,  at  wall 


points  on  edges  1  or  3,  the  quantities  5 


x? 


,,  ,  r  ,  are  required, 

■yc’  ’z?*  M 


Similarly,  n  ,  n  ,  n  ,  must  be  calculated  at  wall  points  on  edges  2  and 

y^> 

4.  All  of  these  quantities  can  be  related  to  the  derivatives  of  T-^  and  I2 
given  in  Eqs.  (3)  and  (4).  The  following  expressions  result: 


=  (Tnsx  "  sn  Tx)/J  '>  =  (-x^sx  +  s^x)/j;  sx  - 


;  =  (t  s  -  s  t  J/j 

y  n  y  n  y' 


ny  *  <-Vy+  Vy)/j;  5y 


Ez  =  <Vc  '  Vs)/;i  '  nz  =  <Vc 


’  r’z  vz 


(A-la) 


K 

'  K 

' 

'  X? 

'y? 
1  z s 


(t  S  +  T  S 
v  n?  X  n  X? 


Vx  -  VxC)/j  -  (5xL)/j 


Vyc)/j  '  (5yV/j 


=  (t  S+tS  -St 

nc  y  n  K  nc  y 

=  <t«T5  +  Vcc  "  VSC  "  V?C)Ai  '  (?zV/j 

=  (t«5x  -  V*  +  scr,Tx  +  W/j  '  ("xV/J 

=  (-T«sy  ■  'csyc  +  Vy  +  Vyc)/J  '  (Vc)/j 

=  (t r  s  +  t_S_  -  S_  T  -  S  t  )  / j  -  (nJJ/j 


(A- lb) 


54  4 


5  44 


54  4 


5  44 


where:  j  =  s_t  -  s  t_ 

5  n  n  5 

j  =  S  T  +  S_T  -s  T_  -  S  T*. 

J4  54  n  5  n4  n4  5  n  54 

sx  =  cos(t  ) ,  sy  =  sin(r),  sz  =  0 

t  =  cos(r)/s,  t  =  -sin(t)/s,  t  =  0 

y  x  z 

ux  =  ey  =  0.  oz  =  1 

s  =  -sin(T)T  ,  s  =  cos(t)t  ,  s  =  0 

x?  4  K  4  zr 

T  =  -[cos(t)s  /s  +  sin(T)Tr]/s 

t  =  [sin(T)s  /s  +  cos(t)t  ]/s  t  =  0 
y4  4  4  zr 


cylindrical  coordinates 


A-l 
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3  l’  sy  3  sz  =  0 


=  1  -  T, 


Tz  =  0 


cartesian  coordinates 


u  =l,u  =u  =0 
z  x  y 

S  =  S  =T  =  T  =0 

x?  ye  ye  xe 


Except  at  points  which  link  adjacent  zones,  the  derivatives  of  t,  s,  with 
respect  to  £,  n,  e»  are  evaluated  analytically.  The  expressions  for  these 
quantities  involve  corner  values  of  these  parameters,  which  can  be  computed 
using  the  formulas  given  below,  for  corner  i: 


where : 


i? 


(q  +  q  0  )/(l  -  q  0  ) 

V  V  ~T  S 


3  ♦  es%>/<1  -  Vs1 


S-  ={q  +q0+0q+s(q0+q0  )}/(l  -  q  0  ) 

mt?  v  veMx  C  S  mt  se';  K  s; 

t  -  =  f0  +  q  0„  +  q  0  +  t  (q  0  +  q  0  )}/(l  -  q  0  ) 

i 1  v?  Mve  s  \  e  'Mxe  s  mt  se'i  '  Hx  s' 


sie  3  (sk  ‘  5j)f'U) 
s i cc  3  (skC  '  sjc)f'({) 


(A-2) 


Tin 


Ti„t  3  (Tmc  -  V9'(n) 


Sin 


9t  (tih  -  ■'PS'(ii) 


Sine  ^"9t  ~  r it)  +  ^ 9ttti 5  +  9T\i^9  ^ n ^ 

T1{  -  es(rk  - 

"  [es(rkc  -  rjC)  +  <rk  -  rj)(0ssr,c  *  ssv):1,"(?> 


i  =  1  => 

q  = 

b. 

=3- 

tl 

© 

J  =  1, 

k  =  2, 

II 

4^ 

m  =  1 

i  =  2  => 

q  = 

c. 

0  =  T, 

j  =  1, 

k  =  2, 

1-3, 

m  =  2 

i  =  3  => 

q  = 

c 

0  =  0, 

j  =  4, 

k  =  3 

cr> 

n 

m  =  2 

II 

-P* 

II 

V 

q  = 

b, 

II 

o 

j  =  4, 

k  =  3, 

£  =  4, 

m  =  1 

In  the  above  equations,  q,  0  and  their  derivatives  are  evaluated  at  the  corner 
under  consideration. 


Expressions  for  the  first  derivatives  of  s,t  with  respect  to  n,  ?  are 
as  follows: 


A-2 
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T-(5.n,c)  =  f'U)  tas[s(5,0,c),v](s3  -  s4)g(n)  +  *s[sU  ,1  ,c)  ,v]g(n)} 
xn(C  n,c)  =  g'(n)  {^[s(?  ,1 ,4) ,v]  -  o[s(c,0  c)  v]} 

T5(5,n.c)  =  g(n)  {osCs(5,0,c)  v]  [s4  ?U)  +  s3  fU)3  +  av[s(?  ,0  .^)  ,vj} 

+  g(n){»s[s(?,l,c)  v]  [si  T(?)  +  s2  f(?)]  +  'F  [s (5 ,1  4hv-l} 

c  c 

sr(5  n  ?)  =  f‘U)  {c[t(?  ,1 ,?)  ,v]  -  b[t(5,0  4)  v]} 

sn(5.n  ?)  =  y1 (n){bT[T(?,0,c) ,v][T1-T4]f(5)+(cTLT(5,0,?) 5v][T2-T3Jf(5)} 

srU  n  ^4)  =  f(c){b  [tU,0  4)  ^v][t4  g(n)+xi  y(n)]+bv[xU,0,4)  ,v]} 

L,  l  -r^ 

+  f (c ) { ct[t( c  ,1 ,?) , v][to  g(n)  +  to  g(n)]  +  c  [t(c,l,?),v]} 

e  4 

where: 


f‘(5)  =  g*(e)  =  ,  T(?)  =  i-f(c).  g(5)  =  i-y(5) 

On  edges  1  and  3,  required  second  derivatives  ot  s,t  with  respect  to  5,n,4> 
are  computed  from: 

t  =  t  g(n)  +  tk  g(n) 

44  44 

V  ‘  <V  -  Tt>9'(n) 

4>  S 


x rr  =  x  g(n)  +  xk  g(n) 
“  54  54 


s  =  ( q  t  +  2q  )  t  +  q  +  q  x 


(A-4) 


nc 


(q  t  +  q  )xr  +  qT^ 


where: 


Srr  =  {cT[<t>2  9(h)  +  <t>3  9(h)]  bx^i  +<l>4  g(n)J  +  (cv_  bv)}t;U) 
5  4  44 

on  edge  1,  then  z  =  4,  k  =  1,  q  =  b 
on  edge  3,  then  £  =  3,  k  =  2,  q=c 


Here  q  and  its  derivatives  are  evaluated  at  the  point  under  considerations. 
Similarly,  required  second  derivatives  on  edges  2  or  4  are: 


\c  ’  sn  f({)  +  sk  fU> 

44  44 


SM  =  (sk  -  h  >f'<0 


(fl-b) 


A-3 
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where: 

Here  0 


S„C  =  S«  >({)  *  Sk  f<«> 
n?  nc 

T«  s  (6SSS<;  *  esv)s5  +  ess55 


V  =  t'1's[s2f(5)  +  "  CTs(s3f^)  +  s4f(^  +  ‘ 


(A-6) 


:  =  (0  s  +  20  )s  +  0  +  0  s 

CC  ss  ?  Sv'  ?  vv  s 


on  edge  2  use  l  =  1,  k  =  2  0=Y 
on  edge  4  use  £  =  4 ,  k=3,  0=ct 

and  its  derivatives  are  evaluated  at  the  point  under  consideration. 
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